As an important prototype to study the reaction involved nitrogen atom, N(²D) + H₂ reaction received considerable attention both experimentally^[1–10] and theoretically,^[11–56] because it plays an important role in atmospheric chemistry, combustion of nitrogen-containing fuels, and explosion processes.

From the experimental side, the nascent vibrational population ratios, the vibrational and rotational state distribution of the N(²D) + H₂ reaction and its isotopic variants as well as the rate constant of the N(²D) + HD isotopic variant were reported by Umemoto et al.^[3–6,8] with improved experimental resolution. In order to produce the N(²D) atom, an intense laser pulse at 275.2 nm was employed to irradiate a mixture of NO and H₂(D₂) and the laser-induced fluorescence was used to detect the products of the N(²D) + H₂ reaction and its isotopic variants. Employing a crossed molecular-beam experiment, Alagia et al.^[7,9] detected the N(²D) + D₂ reaction by mass spectrometric and a forward-backward symmetry differential cross section (DCS) was observed at collision energy 165 meV and 220 meV.

Numerous theoretical calculations have been carried out on the N(²D) + H₂ system. With respect to the potential energy surface aspect, using multi-reference configuration interaction (MRCI) plus Davidson correction and aug-cc-pVTZ basis set, Pederson et al.^[14] constructed a PES for N(²D) + H₂ reaction. Ho et al.^[23] improved the PES which was reported by Pederson et al.^[14] However, these PESs are not accurate enough, the rate constant based on these PESs are always overestimated. On the other hand, the vibrational levels of NH₂ are strongly affected by the Renner–Teller(RT) coupling, which is not included in previous theoretical calculation. So, many groups^{[25,38,39,41,46,47]} have done lots of work for further improving the PES over the next few decades.

For the dynamic aspect, reaction probabilities, integral cross section (ICS), DCS, rate constant, and isotope effect, etc. have been investigated by quasi-classical trajectory (QCT),^{[7,9,11–14,20,21,24,26,33,34,45,48,50,51,55]} quantum dynamics calculation,^{[12,15,20,22,27–37,39,42–45,49,52,56]} and statistical quantum mechanical method.^{[18,24,26,40]} However, there are few works investigating the isotope effects of the N(²D) + HD reaction at the state-to-state level of theory. In order to obtain more accurate information and an accurate physical understanding of the title reaction, the time-dependent quantum wave packet (TDWP) method is employed in the state-to-state dynamic calculation based on the PES of Li et al.^[46] in the collision energy range from 0 eV to 1.0 eV. This article is organized as follows. In Section 2, the theoretical method is briefly introduced. In Section 3, the results and discussion are presented, and the conclusions are presented in Section 4.

The Hamiltonian of the N(²D) + HD reaction in body-fixed (BF) reactant Jacobi coordinates for a given total angular momentum J can be written as

In a wave packet calculation, the initial condition of the wave packet has to be set up before its propagation. The initial wave packet in space-fixed (SF) representation (v₀, j₀, l₀) can be constructed by using the Gaussian wave packet in the R direction and it can be written as

To propagate the initial wave packet in the BF frame, one should transform |JMj₀l₀ɛ〉 to its BF representation counterpart as

During the propagation, the fast Fourier transform method is adopted to act as the radical kinetic energy operator onto the wave packet combination with an L-shaped grid.^[58] The generalized discrete variable representation (DVR) is used to evaluate the action of the potential energy operator, in which the wave packet is converted from the angular finite basis representation (FBR) to a grid representation. To avoid the wave packet reflecting back from the boundaries, a damping function is employed with the same form as that of Ref. [59].

The obtained BF scattering time-independent wave function from an initial state (v₀, j₀, l₀) is obtained by an orthogonal transformation matrix as

Finally, the scattering matrix in the SF frame is transformed into the helicity representation by the standard transformation

Li et al.^[46] constructed a novel global PES of N + H₂ reaction by using the multi-reference configuration interaction (MRCI) approach. The major features of PES are the bond stretching N–H–H and H–N–H configuration, and these features are displayed in Fig. 1. As shown in Fig. 1, for a collinear N–H–H configuration, the C_∞v saddle point is found at about R_HH = 1.532a₀, R_NH = 4.414a₀ and the saddle point with a potential barrier height of 0.222 eV (5.12 kcal/mol). For the linear H–N–H stretch model, a D_∞h linear saddle point is located at R_NH = 1.8716a₀ and a potential well 4.128 eV (95.19 kcal/mol) below the energy of the N(²D)–H₂ asymptote. The adequacy of the PES for dynamic calculations has been demonstrated in Ref. [51].

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. Contour plot of the PES for bond stretching N–H–H (a) and H–N–H (b) linear configuration in internal coordinates (energies are given in unit eV).

To obtain the optimal numerical parameters, many convergence tests have been carried out for J = 0 as summarized in Table 1. Meanwhile the state-to-state calculations for all J up to 50 are conducted to obtain the convergent ICSs and DCSs of the title reaction for collision energies up to 0.5 eV.

Table 1.

Table 1.

Table 1. Numerical parameters used in the dynamic calculations (all parameters are in atomic unit unless otherwise stated). . Grid range and size: R ∈ [0.1,25] r ∈ [0.1,22.5] Angular basis j: [0,95] Initial wave packet: E0 = 0.5 eV δ = 0.18 R0 = 15.0 Matching plane: Rv0 = 16 Absorbing potential: n = n′ = 2 Ra = 21.5 Rb = 24 ra = 20.4 rb = 22.4 Total propagation time: 2000 iterations with Δt = 15.0

Table 1. Numerical parameters used in the dynamic calculations (all parameters are in atomic unit unless otherwise stated). .

3.1. Reaction probability

The total and vibrational resolved reaction probabilities of the two channels of N(²D) + HD reaction for angular momentum J = 0 are shown in Fig. 2. As shown in Fig. 2, the reaction probabilities of the two channels have an intrinsic threshold at E_c ≈ 0.02 eV. This threshold is smaller than the barrier (0.0685 eV), which is located on the minimum energy reactive path. This may imply that the tunneling effect is important in this reaction. Despite the dominance of the deep well, the peaks of total and vibrational resolved reaction probabilities are quite broad. This may indicate that the transition states are short lived for the large exothermicity. For the NH + D channel, the reaction probabilities rise sharply near the threshold energy and increase gradually as the collision energy further increases. Whereas for the ND + H channel the situation is opposite, the reaction probabilities decrease slightly as the collision energy increases. Such phenomena may be aroused by the different zero point energy of NH and ND molecules. The vibrational resolved reaction probabilities of the two channels are also collected in Fig. 2 as a function of the collision energy. All the vibrational state population decreases monotonically with increasing vibrational quantum number in the collision energy from 0 eV to 1.0 eV. From the Fig. 2, the vibrational state resolved reaction probabilities are helpful to determine the contribution of each peak on the total reaction probability.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. Total and vibrational state resolved reaction probabilities of NH + D and ND + H channels are collected in the left and right panel, respectively.

The rotationally resolved reaction probabilities for the two channels are listed in Fig. 3 as a function of collision energy. As shown in Fig. 3, there are many resonances which were dominated by the deep potential well in the entrance channel. For the NH + D channel, it seems irregular for different product rotational states. Whereas for the ND + H channel, the product rotational state j′ from 0 to 3 has similar peaks. Besides, the amplitude of the NH + D channel is larger than the ND + H channel, which may be aroused by the different energy gap between NH and ND molecules.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. The rotational state resolve reaction probabilities of NH + D and ND + H channels for the initial state (v = 0, j = 0) are collected in the left and right panel, respectively.

3.2. Integral cross section

The total and vibrationally resolved ICSs of the two channels of the N(²D) + HD reaction are collected in Fig. 4 as a function of collision energy. As shown in Fig. 4, the total ICS of the NH + D channel is monotonically increased with the collision energy. The ICS of ND + H is also increasing in the collision range from 0 eV to 0.3 eV and it goes into a plateau as the collision energy further increases. Summing populations of rotational states, the vibrationally resolved ICSs of the two channels of the title reaction are also shown in Fig. 4. The product vibrational state population decreases with increasing the vibrational quantum number. This may imply that more rotational states survived in low-lying vibrational states. There are also some mild resonances which can be attributed to the deep well in the entrance channel and large exothermal feature.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. The total and vibrational state resolved integral cross sections of NH + D and ND + H channels are shown in panels (a) and (b), respectively.

The rovibrationally resolved ICSs of the two channels are collected in Figs. 5 and 6 for several selected collision energies. All the rotational distributions are inverted with a peak near the highest rotational states. As shown in Figs. 5 and 6, the number of product vibrational states is almost not changed as the collision energy increases, which may imply that the total energy is mainly transferred to the product translation energy. Comparing Figs. 5 and 6, the ND + H channel has more vibrational and rotational state numbers.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. The rotational state distribution of the NH + D channel at collision energies of 0.038, 0.101, 0.200, and 0.470 eV.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. The rotational state distribution of the ND + H channel at collision energies of 0.038, 0.101, 0.200, and 0.470 eV.

3.3. Differential cross section

The differential cross sections of the two channels are displayed in panels (a) and (b) of Fig. 7 for some selected collision energies. Both angular distributions are peaked at two extreme angles (0° and 180°) at all collision energies. All the figures possess forward–backward symmetry, which may indicated that the insertion mechanism is dominated in the reaction and the statistical model can work well for the title reaction.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. The differential cross sections of NH + D and ND + H channels for some selected collision energies.

3.4. Thermal rate constant

The initial state specific thermal rate constants of the N(²D) + HD reaction are plotted in Fig. 8 along with available experimental data.^[6] The 1²A″ state is just one of the five resulting from the degeneracy of the ²D state of the N atom. So, a statistical factor of 1/5 was considered in the rate constant calculation. The thermal rate constant contributed by rotational excited states need to be thermally average included. However, as the previous theoretical reports^[27,28,30] the contribution by the rotational excited state is very small. Although only the initial state with j₀ = 0 was considered, it may be closer to the thermal rate constant.

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. Initial state specific reaction rate constant for the two channels of N + HD (v = 0, j = 0) as well as the experimental data obtained from Ref. [6].

For the two output channels of the N(²D) + HD reaction, the experimental data are very scare. The sum of the theoretical rate constant of the two output channels is 1.88 × 10⁻¹² cm³ · s⁻¹, which is in excellent agreement with the experimental data^[6] (1.83±0.12 × 10⁻¹² cm³ · s⁻¹) at the temperature T = 300 K. It may imply that the rate constants of the two output channels of the N(²D) + HD reaction are both accurate and reliable.

The N(²D) + HD reaction with initial state (v = 0, j = 0) was calculated by the time-dependent wave packet method in the collision energies from 0 eV to 1.0 eV. The rovibrationally reaction probabilities are investigated for total angular momentum J = 0. Some mild resonances have been found, which can be attributed to the deep well in the entrance channel and large exothermicity. The DCSs have forward–backward symmetry, which may indicate that the insertion mechanism is dominated in the reaction. There are no experimental data of the rate constant for the two channels of the title reaction. However, the sum of the theoretical rate constant of the two output channels are in good agreement with the experimental data. This indicated that the rate constants of the two output channels of the N(²D) + HD reaction may be both accurate and reliable.